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Abstract 

During  directional  solidification  of  a binary  alloy  at  constant  velocity,  thermosolutal 
convection  may  occur  due  to  the  temperature  and  solute  gradients  associated  with 
the  solidification  process.  For  vertical  growth  in  an  ideal  furnace  (lacking  horizontal 
gradients)  a quiescent  state  is  possible.  For  a range  of  processing  conditions,  the 
thermal  Rayleigh  number  is  sufficiently  small  that  the  stabilizing  role  of  the  thermal 
field  during  growth  vertically  upwards  may  be  neglected,  and  only  solutal  convection 
need  be  considered.  The  effect  of  a time-periodic  vertical  gravitational  acceleration  (or 
equivalently  vibration)  on  the  onset  of  solutal  convection  is  calculated  based  on  linear 
stability  using  Floquet  theory.  We  find  that  a stable  base  state  can  be  destabilized 
due  to  modulation,  while  an  unstable  state  can  be  stabilized.  The  flow  and  solute 
disturbance  fields  show  both  synchronous  and  subharmonic  temporal  response  to  the 
driving  sinusoidal  modulation. 
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1.  Introduction 


During  solidification  of  a binary  alloy,  temperature  and  solute  gradients  are  inherently 
present,  and  cause  a density  gradient  in  the  melt.  The  action  of  a gravitational  field  on 
this  density  gradient  can  give  rise  to  fluid  flow  in  the  melt  with  consequent  redistribution  of 
solute,  which  may  result  in  an  undesirable  solute  distribution.  The  possibility  of  processing 
materials  in  space  is  attractive,  since  the  low  level  background  gravitational  acceleration 
(10-6  of  the  earth’s  gravitational  acceleration,  ge ) can  potentially  eliminate  buoyancy-driven 
convection.  However,  owing  to  orbital  maneuvers  and  inherent  mechanical  vibrations,  the 
occurrence  of  time-dependent  local  acceleration  (g-jitter)  may  by  itself  induce  buoyant  con- 
vection. The  local  residual  accelerations  in  space  laboratories  have  been  characterized  [1,  2]. 
On  earth  vibration  may  be  useful  as  a means  for  controlling  convection.  This  idea  has  been 
investigated  recently  for  two  systems  (CsCdCl3  [3]  and  CdTe  [4]),  where  low  frequency  vi- 
brations were  used  in  an  effort  to  control  the  interface  shape  and  position  during  Bridgman 
growth. 

For  vertical  Bridgman  growth,  the  growth  direction  is  aligned  with  the  gravity  vector. 
Ideally,  only  the  vertical  thermal  and  solutal  gradients  from  the  solidification  process  exist 
and  a quiescent  state  in  the  melt  is  possible.  For  growth  of  a pure  material  vertically  upwards, 
the  temperature  will  increase  with  height,  and  for  a normal  melt  which  expands  on  heating 
the  density  will  decrease  with  height.  Hence,  for  growth  vertically  upwards  the  temperature 
field  is  stabilizing  with  respect  to  convective  instabilities.  For  an  alloy  which  rejects  a light 
solute  upon  solidification,  the  solute  gradient  which  develops  above  the  interface  tends  to 
increase  the  density  of  the  melt  with  height.  As  a result  of  this  situation,  either  solutal  or 
thermosolutal  convection  may  occur  depending  on  the  relative  magnitude  of  the  temperature 
and  concentration  gradients. 

The  onset  of  thermosolutal  convection  during  directional  solidification  in  a constant  grav- 
itational field  was  studied  by  Coriell  et  al.  [5],  and  the  onset  of  solutal  convection  alone  by 
Hurle  et  al.  [6].  The  review  article  by  Glicksman  et  al.  [7]  gives  a comprehensive  bibliography 
on  convection  in  directional  solidification,  which  includes  studies  of  the  related  interfacial  in- 
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stability  problem  as  well  as  the  interaction  of  convective  and  interfacial  instabilities.  Results 
from  these  previous  studies  indicate  that  for  dilute  alloys  under  processing  conditions  where 
interfacial  instabilities  should  not  occur,  it  may  be  very  difficult  to  avoid  convection  under 
terrestrial  conditions.  In  order  to  evaluate  the  practicality  of  space-based  processing  using 
the  Bridgman  technique,  the  effect  of  local  time-dependent  residual  accelerations  on  buoyant 
convection  must  be  investigated.  A recent  review  by  Alexander  [8]  discusses  the  effect  of 
residual  accelerations  on  heat  and  mass  transfer  in  low-gravity  materials  experiments.  The 
concern  of  the  present  study  is  whether  sinusoidal  time-dependent  accelerations  can  lead 
to  the  onset  of  convection  in  a quiescent  melt  when  the  accelerations  are  aligned  with  the 
growth  direction. 

The  stability  of  fluid  systems  undergoing  time-periodic  forcing  has  been  reviewed  by 
Davis  [9]  and  Ostrach  [10].  The  bulk  of  the  existing  work  has  concentrated  on  Rayleigh- 
Benard  convection  subject  to  gravity  modulation  or  boundary  temperature  modulation,  and 
Taylor-Couette  flow  with  time-periodic  torsional  oscillation  of  the  cylinders.  Some  early 
work  on  temporally  modulated  Rayleigh-Benard  convection  by  Gershuni  and  Zhukhovitskii 
is  summarized  in  their  book  on  convective  stability  [11].  They  considered  different  aspects 
of  time-periodic  modulation,  including  Floquet  theory  for  the  linear  onset  conditions  and 
finite-difference  calculations  of  the  nonlinear  behavior  for  sinusoidal  modulation.  They  found 
that  both  stabilization  and  destabilization  can  be  obtained  as  a result  of  the  modulation. 
They  also  employed  averaging  techniques  to  study  the  asymptotic  behavior  in  the  case  of 
high  frequency  modulation.  For  Rayleigh-Benard  convection  subject  to  sinusoidal  gravita- 
tional acceleration,  Gresho  and  Sani  [12]  performed  an  extensive  study  of  the  linear  theory 
parameter  space  using  Floquet  analysis  and  a highly  truncated  Galerkin  expansion  for  the 
disturbance  solutions.  In  the  extreme  case  of  a single  expansion  function,  the  analysis  is 
reduced  to  consideration  of  the  stability  behavior  of  the  Mathieu  equation.  They  found 
that  sinusoidal  gravity  modulation  has  a significant  effect  on  the  stability  limits  of  Rayleigh- 
Benard  convection.  Disturbance  response  having  the  same  frequency  as  the  modulation 
(synchronous  behavior)  is  found  for  low  to  moderate  modulation  frequencies,  while  response 
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with  half  the  modulation  frequency  (subharmonic  behavior)  is  obtained  for  high  frequency 
modulation.  In  addition  to  linear  theory,  they  studied  the  nonlinear  behavior  resulting  from 
sinusoidal  modulation  using  a low-order  Galerkin  expansion. 

More  recently,  Wadih  and  Roux  [13]  have  considered  the  effect  of  gravity  modulation  on 
an  infinite  vertical  cylinder  of  fluid  subject  to  a constant  vertical  temperature  gradient.  A 
three-dimensional  simulation  of  Rayleigh-Benard  convection  with  gravitational  modulation 
has  been  performed  by  Biringen  and  Peltier  [14].  Nonlinear  convection  was  simulated  re- 
sulting from  both  sinusoidal  and  random  modulation  with  and  without  a one-<7  background. 
They  obtained  synchronous  and  subharmonic  response  resulting  from  sinusoidal  modulation 
in  agreement  with  Gresho  and  Sani,  and  showed  how  the  finite-amplitude  solutions  transport 
heat  under  various  modulation  conditions. 

In  the  following  section,  the  simplified  model  of  solutal  convection  in  directional  solid- 
ification is  presented.  The  emphasis  of  the  present  work  is  the  conditions  for  the  onset  of 
convection  from  the  quiescent  state  based  on  linear  stability  theory.  Because  of  the  periodic 
time-dependent  forcing,  Floquet  theory  is  used  in  the  stability  analysis.  The  solution  for  the 
onset  conditions  requires  numerical  solutions  of  the  resulting  differential  eigenvalue  problem. 
Two  distinct  numerical  approaches  are  used  to  obtain  the  linear  stability  solutions.  Results 
for  a range  of  parameters  corresponding  to  different  alloy  systems  are  presented.  Some  brief 
results  from  nonlinear  calculations  are  included  to  elucidate  the  hnear  results  and  provide 
an  additional  verification  of  the  hnear  theory  results. 

2.  Governing  equations 

The  model  for  solutally-driven  buoyant  convection  in  directional  solidification  is  shown  in 
Fig.  1.  The  exponential  solute  profile  in  the  melt  results  from  the  rejection  of  solute  at  the 
interface.  Depending  on  the  densities  of  the  alloy  constituents  and  the  growth  direction,  a 
statically  unstable  density  stratification  may  result  from  the  concentration  variation.  For 
steady-state  conditions,  the  solid  composition  is  equal  to  the  far-field  concentration  in  the 
liquid  c^o.  The  solute  jump  at  the  interface  is  characterized  by  the  equilibrium  segregation 
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coefficient  k,  which  is  the  ratio  of  the  solid  and  liquid  concentrations  at  the  interface;  for 
k less  than  one,  solute  is  rejected  at  the  interface.  For  growth  vertically  upwards,  rejection 
of  a light  solute  produces  the  unstable  stratification.  For  the  present  study,  the  effects  of 
the  thermal  field  and  deformation  of  the  interface  will  not  be  included,  but  are  part  of  an 
ongoing  investigation  of  this  problem.  We  also  assume  that  the  densities  of  the  crystal  and 
the  melt  are  equal,  thus  neglecting  the  possibility  of  convection  due  to  a density  change. 

In  order  to  investigate  the  onset  and  development  of  buoyancy-driven  solutal  convection 
during  directional  solidification,  the  governing  equations  are  written  in  terms  of  a Cartesian 
coordinate  system  (x,  y,z)  which  is  attached  to  the  planar  crystal-melt  interface  moving 
at  constant  growth  velocity  V in  the  positive  2 direction.  The  fluid  velocity  in  the  melt 
is  described  by  the  vector  field  u = (u,v,w).  This  velocity  is  measured  in  the  laboratory 
frame  in  which  the  crystal  is  at  rest,  so  that  in  the  undisturbed  state  u = 0.  Assum- 
ing an  incompressible  Newtonian  fluid,  the  basic  governing  equations  for  the  problem  are 
the  conservation  of  mass  equation,  the  Navier-Stokes  equations  written  in  vector  form,  the 
convection-diffusion  equation  for  the  solute  concentration  c,  and  a linear  equation  of  state 
expressing  the  dependence  of  the  melt  density  p on  solute  concentration, 


V • u = 0 

dU  r. 

— + (u*  • V)u  = -Vp/po  + uV2u  + gp/ po 

~ + (u*  ■ V)c  = DV2c 
dt  v ' 

p = p0[l  - (3{c-c0)\ 


(la) 

(lb) 

(lc) 

(ld) 


where  the  velocity  u*  = u — (0,  0,  — V)  results  from  fixing  the  reference  frame  to  the  moving 
solid-liquid  interface.  Here  p is  the  fluid  pressure,  v is  the  kinematic  viscosity,  D is  the  solute 
diffusion  coefficient,  and  (3  is  the  solutal  expansion  coefficient.  The  Boussinesq  approximation 
is  employed  which  assumes  that  density  variation  is  important  only  in  the  body  force  term 
in  the  Navier-Stokes  equations.  The  remaining  fluid  parameters  {v,  D,  and  (3)  are  assumed 
constant.  The  subscript  zero  refers  to  reference  values  of  the  parameters. 
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For  a gravitational  field  aligned  with  the  growth  direction,  a motionless  solution  exists. 
We  consider  the  linear  stability  of  this  motionless  base  state  subject  to  time-periodic  gravita- 
tional acceleration.  The  model  for  solutal  convection  during  directional  solidification  assumes 
growth  with  a rigid  planar  interface  in  the  presence  of  a vertical  gravitational  field  that  varies 
sinusoidally  in  time  according  to  g(t)  = go  + g\  cos  (Of),  where  Q is  the  dimensional  forcing 
frequency.  The  liquid  is  assumed  to  be  semi-infinite  in  the  vertical  (2)  direction,  and  laterally 
unbounded  so  that  solutions  that  are  periodic  in  the  horizontal  directions  are  possible.  The 
steady-state  base  solution  consists  of  no  fluid  motion  and  the  exponential  solute  profile  ob- 
tained from  solving  the  steady,  one-dimensional  diffusion  equation  in  the  semi-infinite  liquid 
region  subject  to  solute  conservation  at  the  growing  interface,  given  by 


D^-  = -V(l-k)c, 


(2) 


and  subject  to  the  far-field  concentration  value  Coo. 

For  the  linear  stability  analysis  of  the  base  state,  the  flow  field  variables  are  written  as 
the  superposition  of  the  base  state  component  and  a perturbation.  The  perturbed  quantities 
are  Fourier  analyzed  in  the  lateral  directions,  so  that  the  flow  variables  are  written  as 


/ f ,V  \ 

u{x,  y,z,  t) 

0 1 

u(z,t ) ^ 

v(x,y,z,t) 

0 

v(z,t), 

w(x,y,z,t) 

= 

0 

+ 

w{z,t) 

p(x,y,z,t) 

?'%-) 

p(M) 

^ c{x,y,z,t)  y 

v c(M)  J 

exp^a^x  + iayy ), 


(3) 


where  ax  and  ay  are  the  wavenumbers  in  the  lateral  directions.  The  base  state  components 
are  denoted  by  a zero  superscript  and  the  quantities  u,  v,  etc  are  the  perturbation  amplitudes. 
Governing  equations  for  the  perturbation  quantities  are  obtained  by  substituting  the  above 
quantities  into  the  set  of  equations  (1)  and  linearizing  in  the  perturbation  quantities.  The 
set  of  linearized  equations  for  the  perturbation  quantities  contain  time-dependent  coefficients 
resulting  from  the  assumed  form  of  the  gravitational  acceleration. 

For  the  linear  stability  analysis,  the  problem  can  be  reduced  to  solving  a single  fourth- 
order  equation  for  the  the  perturbed  vertical  velocity  w,  using  the  standard  manipulation  for 
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hydrodynamic  stability  analyses  to  remove  the  pressure  [15],  and  the  second-order  equation 
for  the  perturbed  solute  concentration  c.  The  variables  are  made  dimensionless  by  introduc- 
ing the  length  scale  D /V  and  the  time  scale  D/V2,  and  by  measuring  the  flow  velocity  in 
units  of  the  solidification  velocity  V and  the  solute  field  in  units  of  the  far  field  concentration 
Coo.  The  dimensionless  linearized  equations  for  the  stability  problem  are  then 

— {Lwt  — Lwz ) = L2w  — 4-  Rs^  cos(flt)]c,  (4a) 

Sc 

ct  - cz  + wc[0)  = Lc , (4b) 

Here,  subscripts  indicate  partial  derivatives,  and  = —[(1  — k)/k]e~z  is  the  gradient  of 
the  unperturbed  solute  field.  We  have  defined  the  operator  L = (d2 /dz2  — a2),  where  a 
is  the  spatial  wavenumber  in  the  horizontal  directions  (a  = yja2  -fi  a2).  Sc  — vfD  is  the 
Schmidt  number,  Rs^  — gofic^D /V)3 / {vD)  is  the  Rayleigh  number  based  on  the  mean  or 
d-c  gravitational  field,  Rs^  = gxfic^D /V)3  / (vD)  is  the  Rayleigh  number  based  on  the  a-c 
part  of  the  gravitational  field,  and  Q.  = QD/V 2 is  the  dimensionless  oscillation  frequency. 
Table  1 summarizes  the  dimensional  and  dimensionless  parameters. 

The  dimensionless  boundary  conditions  at  the  interface,  2 = 0,  are 


w = wz  = 0 

(5a) 

cz  + (1  — k)c  = 0. 

(5b) 

The  perturbations  are  required  to  decay  as  2 — *■  oo.  For  numerical  purposes,  we  assume  that 
the  spatial  domain  extends  from  the  interface  2 = 0 to  a truncated  value  of  the  semi-infinite 
domain  denoted  by  hi.  The  perturbation  quantities  are  set  equal  to  zero  at  the  far-field 
boundary. 

Solutions  to  the  above  set  of  equations  which  have  time-periodic  coefficients  can  be 
obtained  using  the  framework  of  Floquet  theory  [16],  in  which  the  solutions  are  represented 
as  the  product  of  a temporally  periodic  function  and  an  exponential  function  of  time.  Two 
different  numerical  implementations  of  Floquet  theory  were  employed  to  solve  the  stability 
problem  formulated  above  and  are  discussed  briefly  below. 
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The  first  approach  consists  of  representing  the  periodic  component  of  the  solutions  by  a 
truncated  complex  Fourier  series  in  time.  The  solutions  for  the  perturbation  quantities  w 
and  c are  represented  by  the  product  of  a periodic  Fourier  series  and  an  exponential  term 
with  complex  growth  rate  cr, 

w(z,t)  = e"‘  £ wm{z)e'm!u, 

| m | < A/ 

c(z,t)  = e"  £ cm(z)e'mn‘. 

|m|<Af 

Substitution  of  the  series  expansion  into  the  equations  and  boundary  conditions  above  yields 
a set  of  12 M + 6 coupled  two-point  boundary  value  problems  in  the  spatial  variable  2 for 
the  complex  Fourier  coefficients  (inm,cm)  of  the  periodic  solution  components  . The  re- 
sulting set  of  coupled  ordinary  differential  equations  subject  to  the  homogeneous  boundary 
conditions  yields  an  eigenvalue  problem  that  is  solved  in  a manner  similar  to  the  one  de- 
scribed in  [5].  The  homogeneous  eigenvalue  problem,  is  converted  into  an  inhomogeneous, 
nonsingular  problem  using  the  approach  suggested  by  Keller  [17].  The  coupled  set  of  linear 
two-point  boundary  value  problems  is  solved  using  the  computer  code  SUPORT  [18],  which 
uses  superposition  of  numerically  integrated  solutions  with  an  orthonormalization  procedure 
to  maintain  the  linear  independence  of  the  solution  set.  A high-order  Adams-type  method 
is  used  for  the  numerical  integration  of  the  spatial  dependence  in  the  SUPORT  code.  Here, 
the  temporal  resolution  depends  on  the  number  of  coupled  Fourier  modes  (A/)  used  to  ap- 
proximate the  periodic  part  of  the  solution.  The  value  of  M depends  most  strongly  on  the 
modulation  frequency,  with  lower  frequencies  requiring  more  temporal  modes  to  obtain  a 
given  level  of  accuracy. 

The  set  of  equations  contains  the  complex  growth  rate  cr  as  a parameter.  The  base 
state  subject  to  periodic  forcing  is  linearly  stable  for  a given  set  of  parameters  if,  for  all 
disturbances,  crr  < 0,  where  crr  is  the  real  part  of  cr.  In  the  calculations,  setting  crr  = 0 
allows  for  the  determination  of  marginal  values  of  the  modulation  amplitude,  Rs^,  and  the 
imaginary  part  of  the  growth  rate,  cr,,  for  fixed  values  of  the  remaining  parameters.  For  the 
results  presented,  only  neutral  disturbances  with  cr,  = 0 or  Q/2  are  obtained.  In  this  solution 
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approach,  the  number  of  complex  valued  unknowns  is  reduced  to  6 M + 6,  if  the  marginal 
values  of  Rs^  are  determined  for  the  specific  values,  crx  = 0 or  fl/2. 

The  second  approach  employed  consists  of  approximating  the  spatial  behavior  of  the 
disturbance  quantities  via  the  pseudospectral  technique  in  the  physical  domain  as  described 
in  [19].  The  approach  corresponds  to  expanding  the  solutions  in  terms  of  truncated  series  of 
Chebyshev  polynomials  Tn(s), 

N 

w(z,t ) = WM  ^n(s), 

n=0 

N 

c{z,t ) = Tn{s), 

n=0 

where  s = (2 z/hi)  — 1.  The  pseudospectral  discretization  requires  that  the  solution  ex- 
pansions satisfy  the  governing  equations  at  specific  collocation  points  for  the  Chebyshev 
polynomials.  When  implemented  in  the  physical  domain  the  unknowns  are  the  solution 
values  at  the  collocation  points.  The  spatial  differential  operators  in  the  governing  partial 
differential  equations  are  replaced  by  discrete  matrix  operators.  As  a result,  the  governing  set 
of  partial  differential  equations  and  boundary  conditions  becomes  a set  of  coupled  ordinary 
differential-algebraic  equations  in  time  for  the  unknown  solution  values  at  the  collocation 
points. 

The  computer  code  DASSL  [20]  is  used  to  solve  the  differential-algebraic  system  for  one 
complete  period  of  the  driving  temporal  modulation.  The  algorithm  uses  backward  differen- 
tiation of  up  to  fifth-order  to  meet  specified  local  error  tolerances.  This  integration  procedure 
is  well-suited  for  stiff  behavior,  which  is  often  exhibited  by  differential-algebraic  systems  [21]. 
In  this  second  solution  approach,  Floquet  analysis  [16]  is  implemented  by  constructing  a fun- 
damental solution  matrix.  The  columns  of  this  matrix  are  linearly  independent  calculated 
solutions  for  the  unknowns  at  the  end  of  one  forcing  period.  The  eigenvalues  of  this  matrix 
are  the  Floquet  multipliers  from  which  the  complex  growth  rate  a is  obtained.  The  utility 
of  the  second  approach  is  that  cr  itself  is  the  eigenvalue,  so  that  its  magnitude  is  obtained 
for  a given  Rs(l\  which  simplifies  the  search  for  marginal  values  and  gives  more  information 
regarding  the  proximity  of  other  modes. 
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In  order  to  check  our  numerical  results,  we  have  used  both  methods  for  selected  cases 
and  obtained  agreement  to  five  significant  figures.  Accurate  results  are  obtained  for  the 
parameter  range  studied  using  from  M = 4 to  8 temporal  Fourier  modes  in  the  first  approach. 
The  required  number  of  modes  increases  as  the  frequency  decreases.  A relative  local  error 
tolerance  of  10-10  is  used  for  spatial  integration.  In  the  second  approach,  the  Chebyshev 
series  is  truncated  with  N = 20.  A convergence  study  indicated  this  number  of  spatial  modes 
is  sufficient  for  the  parameter  range  studied.  For  the  temporal  integration,  a relative  local 
error  tolerance  of  10-6  was  chosen.  These  same  two  solution  methods  were  used  for  a related 
study  of  modulated  Taylor-Couette  flow,  where  additional  details  concerning  accuracy  and 
convergence  are  given  [22]. 

As  an  additional  check,  the  full  two-dimensional  nonlinear  equations  were  solved  using 
the  approach  described  in  [23].  Briefly,  the  flow  field  is  calculated  in  two  dimensions  using 
the  vorticity-streamfunction  approach,  which  eliminates  the  need  to  calculate  the  pressure 
explicitly.  Standard  second-order  finite-difference  techniques  are  used  to  solve  for  the  finite 
amplitude  behavior  of  the  two-dimensional  velocity  and  solute  concentration  fields. 

3.  Results  and  Discussion 

Even  for  the  present  simplified  model  of  directional  solidification,  a significant  number  of 
parameters  appear  in  the  problem.  For  the  initial  calculations,  conditions  relevant  to  a 
microgravity  processing  environment  are  assumed.  With  a background  gravitational  acceler- 
ation of  10~6<?e,  the  steady-state  solutal  Rayleigh  number  Rs^  will  typically  be  very  small. 
It  is  of  interest  to  determine  whether  sinusoidal  acceleration  in  time  can  lead  to  the  onset 
of  solutal  convection.  So  initially,  we  set  Rs^  = 0,  and  look  for  values  of  Rs^^  that  lead 
to  instability  for  given  values  of  the  acceleration  frequency.  Calculations  are  performed  for 
Schmidt  numbers  10  (corresponding  to  semiconductor  type  alloys),  81  and  144;  the  latter 
two  values  correspond  to  the  lead-tin  and  tin-bismuth  systems,  respectively.  The  value  k = 
0.3  is  used  for  the  lead-tin  system  and  k = 0.28  is  used  for  tin-bismuth.  The  extent  of  the 
liquid  region,  hi  is  fixed  at  15  for  all  the  results  presented;  this  value  is  sufficiently  large  that 


-10- 


it  approximates  a semi-infinite  domain. 

In  Fig.  2,  we  show  the  neutral  stability  boundaries  (crr  = 0)  at  a fixed  value  of  the  lateral 
wavenumber,  a = 0.5,  for  Sc  = 10,  k = 0.3,  hi  = 15,  and  Rs^  = 0.  The  value  a = 0.5 
was  chosen,  since  it  is  close  to  the  critical  value  for  the  unmodulated  case.  The  modulation 
amplitude,  Rs^  is  plotted  versus  the  inverse  of  the  modulation  frequency.  Unstable  regions 
are  shaded.  For  the  range  of  frequencies  investigated,  two  unstable  regions  are  obtained  cor- 
responding to  growing  disturbances  exhibiting  either  synchronous  or  subharmonic  temporal 
response.  At  higher  frequencies,  the  the  theory  predicts  unstable  subharmonic  disturbances, 
where  <7,  = Si/ 2.  The  subharmonic  stability  boundary  is  lobe-shaped,  and  approaches  infin- 
ity as  the  modulation  frequency  becomes  large.  The  synchronous  stability  boundary  (where 
<7,  = 0)  lies  to  the  right  of  the  subharmonic  lobe,  but  there  is  a small  overlap  region.  The 
complete  structure  of  the  synchronous  region  was  not  investigated.  The  interest  here  is  in 
the  minimum  modulation  amplitude  that  yields  instability.  It  is  entirely  possible  that  a 
banded  structure  of  alternating  modes  exists  above  the  lowest  synchronous  boundary;  this 
is  a known  characteristic  of  some  simpler  modulated  stability  problems  [11]. 

It  is  interesting  to  compare  the  magnitude  of  the  modulation  amplitude  required  to  cause 
instability  with  the  critical  solutal  Rayleigh  number  for  the  unmodulated  case.  From  the 
results  of  Hurle  et  al.  [6]  for  Schmidt  number  10,  the  critical  solutal  Rayleigh  number  is 
4.43.  Thus,  for  a modulation  amplitude  of  approximately  fifteen  times  the  steady  critical 
value,  instability  resulting  from  sinusoidal  modulation  occurs  for  a range  of  frequencies.  Fig. 
3 shows  the  same  type  of  stability  diagram  for  Sc  = 81  with  all  remaining  parameters  held 
fixed.  The  diagram  is  qualitatively  similar,  but  both  the  horizontal  and  vertical  scales  have 
changed  substantially.  The  minimum  of  the  subharmonic  branch  is  larger  than  the  Schmidt 
number  10  value,  and  its  width  is  smaller.  The  dimensional  modulation  period  is  2tt  D / {V2Sl)\ 
for  diffusion  coefficients  on  the  order  of  10 ~hcm2  / s and  crystal  growth  velocities  on  the  order 
on  10 _4cm/s,  a value  Si  = 10  corresponds  to  a period  of  600  seconds.  Our  primary  interest 
is  in  shorter  time  periods,  so  we  have  not  extended  the  calculations  to  smaller  values  of  Si. 

Results  representing  the  tin-bismuth  system,  for  which  space  experiments  are  planned 


-11- 


[24],  are  shown  in  Fig.  4.  For  this  alloy,  Sc  = 144  and  k = 0.28.  The  qualitative  behavior  is 
the  same  as  Figs.  2 and  3.  For  this  higher  Schmidt  number,  the  trend  observed  for  Schmidt 
number  10  and  81  continues,  with  the  minimum  of  the  subharmonic  mode  increasing  and  the 
width  of  the  lobe  decreasing.  The  reduction  of  k from  0.3  to  0.28  in  this  set  of  calculations 
has  some  effect  on  the  magnitude  of  the  modulation  amplitude,  e.g.,  with  fi  = 35,  Rs W = 
570.1  and  517.4  for  k = 0.3  and  0.28,  respectively. 

All  the  results  presented  so  far  are  for  fixed  lateral  wavenumber.  To  actually  determine 
the  stability  of  the  system,  one  must  find  the  minimum  value  of  Rs  W for  all  values  of  a.  We 
have  not  done  this  in  general  due  to  the  extensive  number  of  calculations  required.  However, 
for  selected  frequencies,  we  have  found  the  minimum  as  a function  of  a,  and  these  results  are 
shown  as  the  solid  dots  in  Fig.  4.  There  is  only  a small  difference  between  the  a = 0.5  values 
and  the  minimum  values,  indicating  that  fixing  a at  the  value  0.5  is  a good  approximation. 

We  have  investigated  in  detail  the  dependence  of  Rs (b  on  a for  a fixed  Q = 32.0.  The 
results  are  shown  in  Fig.  5.  There  is  subharmonic  branch  at  low  values  of  a and  a synchronous 
branch  at  higher  values  of  a,  with  the  two  branches  intersecting  in  the  vicinity  of  a = 0.45.  For 
a = 0.5,  as  RsW  increases  the  system  becomes  unstable  to  a subharmonic  disturbance,  then 
regains  stability  and  finally  becomes  unstable  to  a synchronous  disturbance.  The  behavior 
is  also  evident  from  Fig.  4 at  this  frequency.  However,  when  one  considers  stability  with 
respect  to  all  possible  values  of  a,  the  system  is  unstable  for  Rs^  greater  than  the  minimum 
value  on  the  subharmonic  branch,  which  occurs  at  a about  0.4. 

The  previous  results  were  for  Rs = 0,  so  that  the  system  is  stable  without  modulation. 
We  also  consider  whether  sinusoidal  modulation  can  stablize  an  otherwise  unstable  system. 
As  mentioned  earlier,  the  critical  Rayleigh  number  for  the  unmodulated  system  is  about 
5,  with  a weak  Schmidt  number  dependence  [6],  We  set  Rs^  = 20,  Q = 200,  and  k = 
0.28,  and  calculate  the  stability  of  the  system  for  Schmidt  numbers  of  10  and  144.  The 
modulation  amplitudes,  Rs^\  as  a function  of  wavenumber,  a,  are  shown  in  Figs.  6 and  7. 
For  small  Rs^\  the  system  is  unstable  to  a synchronous  mode.  However,  as  the  modulation 
amplitude  increases,  the  system  becomes  stable  for  a range  of  Rs^ , losing  stability  again  to 
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a subharmonic  mode  at  sufficiently  large  values  of  Rs^ . The  calculations  were  limited  to  a 
greater  than  0.25.  Disturbances  with  a less  than  0.25  correspond  to  wavelengths  greater  than 
25  D/V,  which  is  on  the  order  of  a typical  container  width.  In  addition,  these  disturbance 
wavelengths  are  larger  than  the  vertical  extent  of  the  domain  (hi  = 15)  which  makes  the 
calculations  sensitive  to  the  value  of  hi.  For  sufficiently  large  values  of  a,  the  system  is  stable 
for  zero  modulation.  The  marginal  curve  (Rs1'0)  vs  a)  for  the  unmodulated  case  is  U-shaped; 
for  Rs = 20,  there  is  a range  of  a for  which  the  system  is  unstable.  For  the  cases  with 
modulation  shown  in  Figs.  6 and  7,  the  region  of  instability  intersects  the  axis  Rs*1*  = 0 (no 
modulation)  in  precisely  this  range.  Because  the  calculations  were  restricted  to  a greater 
than  0.25,  only  the  upper  limit  of  this  range  is  shown. 

To  help  illustrate  the  temporal  character  of  the  disturbances,  we  have  used  the  nonlinear 
code  described  previously  to  compute  finite  amplitude  solutions  for  Q = 20  and  three  values 
of  Rs , as  shown  in  Fig.  8.  The  solutions  are  computed  over  nine  periods  of  the  gravitational 
acceleration,  and  the  growth  or  decay  of  the  solution  is  examined  to  compare  with  the  results 
of  the  linear  analysis  shown  in  Fig.  3.  We  have  plotted  the  value  of  the  streamfunction  at 
a fixed  point  in  the  domain  as  a function  of  time  for  the  gravitational  acceleration  g shown 
in  the  lower  subplot.  For  Rs^  = 400,  the  system  is  unstable  with  a period  twice  that  of 
the  modulation  (subharmonic  response),  in  agreement  with  the  prediction  of  linear  theory. 
If  the  solutal  Rayleigh  number  is  increased  to  Rs^  = 500,  the  system  is  stable,  with  a 
slow  decay  in  the  maximum  values  assumed  by  the  streamfunction  over  two  periods  of  the 
modulation.  The  linear  analysis  for  Rs1'1)  = 500  shows  that  the  system  is  linearly  stable, 
with  a slowest  decaying  mode  that  is  subharmonic;  there  is  also  a synchronous  decaying 
mode  with  a decay  rate  that  is  three  times  larger  than  the  subharmonic  mode.  Finally,  for 
Rs1'1'*  = 600  the  system  is  unstable  with  a synchronous  response,  which  is  also  in  agreement 
with  the  prediction  of  the  linear  analysis. 

In  Fig.  9 we  plot  the  linear  eigenfunctions,  w(z , t)  and  c(z,  f),  at  the  onset  of  instability  of 
the  subharmonic  mode  for  D = 20  and  the  parameters  of  Fig.  3.  The  modulation  amplitude, 
Rs™  is  339.6.  The  spatial  coordinate  is  normalized  by  hi  and  time  is  normalized  by  the 
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modulation  period.  Since  this  is  a subharmonic  mode,  its  period  is  twice  the  modulation 
period.  The  decay  of  the  eigenfunctions  as  2 — > 1 shows  that  the  value  of  hi  used  approxi- 
mates the  semi-infinite  region  well.  The  eigenfunctions  for  subharmonic  response  exhibit  the 
symmetry  f(z,t  + T)  = —f(z,t),  where  T is  the  modulation  period. 

In  summary,  in  the  absence  of  mean  gravitational  acceleration  ( Rs = 0)  and  for  large 
frequencies,  instability  occurs  with  subharmonic  temporal  response  to  the  driving  modula- 
tion. The  amplitude  of  modulation,  Rs^l\  necessary  for  instability  increases  with  increasing 
Schmidt  number  and  modulation  frequency.  At  lower  frequencies,  instability  occurs  with 
synchronous  response.  When  the  system  is  unstable  in  the  absence  of  modulation,  there  is 
a range  of  modulation  amplitudes  for  which  the  system  is  stable. 

In  the  absence  of  modulation  the  critical  Rayleigh  number  for  the  onset  of  convection 
is  only  weakly  dependent  on  the  Schmidt  number  for  large  Schmidt  numbers  [6].  However, 
for  the  case  of  pure  modulation  (Rs^  = 0),  it  is  clear  from  Figs.  2-4  that  Rs1'1'1  depends 
strongly  on  the  Schmidt  number. 
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Table  I 


Definition  of  Symbols 


Dimensional  Parameters 


liquid  diffusion  coefficient 

D 

cm2/s 

solidification  velocity 

V 

cm/s 

bulk  liquid  concentration 

Coo 

wt% 

kinematic  viscosity 

V 

cm2/s 

solutal  expansion  coefficient 

P 

wt%  — 1 

vertical  domain  height 

H 

cm 

gravitational  acceleration 

g(t)  = go  + gi  cos  {Sit) 

cm/s2 

Dimensionless  Parameters 
distribution  coefficient  k 

Schmidt  number  Sc  = v / D 

unmodulated  solutal  Rayleigh  number  Rs^  = gofic^D /V)3/(uD) 
modulated  solutal  Rayleigh  number  Rs^  = gl/3c00(D/V)3/(vD) 
oscillation  frequency  Q = SlD /V2 

vertical  domain  height  hi  = HV / D 
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Figure  Captions 


Figure  1.  Schematic  diagram  of  directional  solidification  model. 

Figure  2.  The  modulation  amplitude  at  the  onset  of  solutal  convection  as  a function  of  the 
inverse  of  the  frequency  for  sinusoidal  gravitational  acceleration. 

Figure  3.  The  modulation  amplitude  at  the  onset  of  solutal  convection  as  a function  of  the 
inverse  of  the  frequency  for  sinusoidal  gravitational  acceleration. 

Figure  4.  The  modulation  amplitude  at  the  onset  of  solutal  convection  as  a function  of 
the  inverse  of  the  frequency  for  sinusoidal  gravitational  acceleration.  Solid  points  indicate 
minimum  with  respect  to  wavenumber  a. 

Figure  5.  The  modulation  amplitude  at  the  onset  of  solutal  convection  as  a function  of  the 
wavenumber. 

Figure  6.  The  modulation  amplitude  at  the  onset  of  solutal  convection  as  a function  of  the 
wavenumber.  In  the  absence  of  modulation  the  system  is  unstable. 

Figure  7.  The  modulation  amplitude  at  the  onset  of  solutal  convection  as  a function  of  the 
wavenumber.  In  the  absence  of  modulation  the  system  is  unstable. 

Figure  8.  Nonlinear  calculations  showing  the  value  of  the  streamfunction  at  a fixed  point 
in  the  melt  as  a function  of  time.  The  time-periodic  gravitational  acceleration  is  shown  in 
the  lower  figure.  The  system  is  unstable  for  Rs = 400  and  600,  and  is  stable  for  Rs1'1'1  = 
500. 

Figure  9.  Linear  eigenfunctions  w(z,  t)  and  c(z,  t)  plotted  for  two  modulation  periods  (sub- 
harmonic response)  for  Q = 20,  Rs = 339.5,  Rs^  = 0,  Sc  = 81,  and  a = 0.5. 
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thermal  Rayleigh  number  is  sufficiently  small  that  the  stabilizing  role  of  the  thermal 
field  during  growth  vertically  upwards  may  be  neglected,  and  only  solutal  convection 
need  be  considered.  The  effect  of  a time-periodic  gravitational  acceleration  (or 
equivalently  vibration)  on  the  onset  of  solutal  convection  is  calculated  based  on 
linear  stability  using  Floquet  theory.  We  find  that  a stable  base  state  can  be 
destabilized  due  to  modulation,  while  an  unstable  state  can  be  stabilized.  The  flow 
and  solute  disturbance  fields  show  both  synchronous  and  subharmonic  temporal  response 
to  the  driving  sinusoidal  modulation. 
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